Libraries

library(tidyverse)
library(reshape2)
library(sf)
library(spdep)
library(tmap)
#library(OasisR) ## old package (OK for version R.2.1)

Synthetic population - H24 library (section 2.1)

Import data for synthetic population located in ‘night’ cell pop0, in ‘day’ cell pop1 and in ‘evening’ cell pop2

pop_0 <- st_read( "data/population.gpkg",
              layer="Population_0", quiet = TRUE)
pop_1 <- st_read("data/population.gpkg",
                  layer="Population_1", quiet = TRUE)
pop_2 <- st_read("data/population.gpkg",
                  layer="Population_2", quiet = TRUE)

Create new colums : gender (2 groups) ; Age (3 groups); Educ (3 Groups)

pop_0 <- pop_0 %>% 
  rowwise() %>% 
  mutate(
    men =  sum(across(starts_with('pop_1'))),
    women = sum(across(starts_with('pop_2'))),
         
    age1 = sum(across(starts_with(c('pop_1_1', 'pop_2_1')))),
    age2 = sum(across(starts_with(c('pop_1_2', 'pop_2_2')))),
    age3 = sum(across(starts_with(c('pop_1_3', 'pop_2_3')))),
    
    educ1 = sum(across(ends_with('_1'))),
    educ2 = sum(across(ends_with('_2'))),
    educ3 = sum(across(ends_with('_3')))
    )

pop_1 <- pop_1 %>% 
  rowwise() %>% 
  mutate(
    men =  sum(across(starts_with('pop_1'))),
    women = sum(across(starts_with('pop_2'))),
    
    age1 = sum(across(starts_with(c('pop_1_1', 'pop_2_1')))),
    age2 = sum(across(starts_with(c('pop_1_2', 'pop_2_2')))),
    age3 = sum(across(starts_with(c('pop_1_3', 'pop_2_3')))),
    
    educ1 = sum(across(ends_with('_1'))),
    educ2 = sum(across(ends_with('_2'))),
    educ3 = sum(across(ends_with('_3')))
  )

pop_2 <- pop_2 %>% 
  rowwise() %>% 
  mutate(
    men =  sum(across(starts_with('pop_1'))),
    women = sum(across(starts_with('pop_2'))),
    
    age1 = sum(across(starts_with(c('pop_1_1', 'pop_2_1')))),
    age2 = sum(across(starts_with(c('pop_1_2', 'pop_2_2')))),
    age3 = sum(across(starts_with(c('pop_1_3', 'pop_2_3')))),
    
    educ1 = sum(across(ends_with('_1'))),
    educ2 = sum(across(ends_with('_2'))),
    educ3 = sum(across(ends_with('_3')))
  )

Distribution of the synthetic population in ’night’ cells in the Paris Region - Maps (Figure 1)

!! TO DO Clementine

Duncan Dissimilarity Index (Figure 1 and Appendix 3)

Duncan’s segregation index is one-group form of dissimilarity index DI and measures the unevenness of a group distribution compared to the rest of the population. It can be interpreted as the share of the group that would have to move to achieve an even distribution compared to the rest of the population. You should not include a column with total population, because this will be interpreted as a group.

In night cells

dfpop_0 <- pop_0 %>% st_drop_geometry()

pop_0_Duncan_18categ<- as.data.frame(ISDuncan(dfpop_0[,2:19]))
colnames(pop_0_Duncan_18categ)[1] <- "Duncan_night"
rownames(pop_0_Duncan_18categ) <- paste (colnames(dfpop_0[ , 2:19]))

pop_0_Duncan_sex<- as.data.frame(ISDuncan(dfpop_0[,20:21]))
colnames(pop_0_Duncan_sex)[1] <- "Duncan_night"
rownames(pop_0_Duncan_sex) <- paste (colnames(dfpop_0[ , 20:21]))

pop_0_Duncan_age<- as.data.frame(ISDuncan(dfpop_0[,22:24]))
colnames(pop_0_Duncan_age)[1] <- "Duncan_night"
rownames(pop_0_Duncan_age) <- paste (colnames(dfpop_0[ , 22:24]))

pop_0_Duncan_educ<- as.data.frame(ISDuncan(dfpop_0[,25:27]))
colnames(pop_0_Duncan_educ)[1] <- "Duncan_night"
rownames(pop_0_Duncan_educ) <- paste (colnames(dfpop_0[ , 25:27]))

pop_0_all_Duncan <- bind_rows(pop_0_Duncan_18categ, pop_0_Duncan_sex, pop_0_Duncan_age, pop_0_Duncan_educ)

pop_0_all_Duncan_round <- pop_0_all_Duncan %>% mutate(across(starts_with("Duncan"), round, 3))

In day cells

dfpop_1 <- pop_1 %>% st_drop_geometry()

pop_1_Duncan_18categ<- as.data.frame(ISDuncan(dfpop_1[,2:19]))
colnames(pop_1_Duncan_18categ)[1] <- "Duncan_day"
rownames(pop_1_Duncan_18categ) <- paste (colnames(dfpop_1[ , 2:19]))

pop_1_Duncan_sex<- as.data.frame(ISDuncan(dfpop_1[,20:21]))
colnames(pop_1_Duncan_sex)[1] <- "Duncan_day"
rownames(pop_1_Duncan_sex) <- paste (colnames(dfpop_1[ , 20:21]))

pop_1_Duncan_age<- as.data.frame(ISDuncan(dfpop_1[,22:24]))
colnames(pop_1_Duncan_age)[1] <- "Duncan_day"
rownames(pop_1_Duncan_age) <- paste (colnames(dfpop_1[ , 22:24]))

pop_1_Duncan_educ<- as.data.frame(ISDuncan(dfpop_1[,25:27]))
colnames(pop_1_Duncan_educ)[1] <- "Duncan_day"
rownames(pop_1_Duncan_educ) <- paste (colnames(dfpop_1[ , 25:27]))

pop_1_all_Duncan <- bind_rows(pop_1_Duncan_18categ, pop_1_Duncan_sex, pop_1_Duncan_age, pop_1_Duncan_educ)

pop_1_all_Duncan_round <- pop_1_all_Duncan %>% mutate(across(starts_with("Duncan"), round, 3))

In evening cells

dfpop_2 <- pop_2 %>% st_drop_geometry()

pop_2_Duncan_18categ<- as.data.frame(ISDuncan(dfpop_2[,2:19]))
colnames(pop_2_Duncan_18categ)[1] <- "Duncan_evening"
rownames(pop_2_Duncan_18categ) <- paste (colnames(dfpop_2[ , 2:19]))

pop_2_Duncan_sex<- as.data.frame(ISDuncan(dfpop_2[,20:21]))
colnames(pop_2_Duncan_sex)[1] <- "Duncan_evening"
rownames(pop_2_Duncan_sex) <- paste (colnames(dfpop_2[ , 20:21]))

pop_2_Duncan_age<- as.data.frame(ISDuncan(dfpop_2[,22:24]))
colnames(pop_2_Duncan_age)[1] <- "Duncan_evening"
rownames(pop_2_Duncan_age) <- paste (colnames(dfpop_2[ , 22:24]))

pop_2_Duncan_educ<- as.data.frame(ISDuncan(dfpop_2[,25:27]))
colnames(pop_2_Duncan_educ)[1] <- "Duncan_evening"
rownames(pop_2_Duncan_educ) <- paste (colnames(dfpop_2[ , 25:27]))

pop_2_all_Duncan <- bind_rows(pop_2_Duncan_18categ, pop_2_Duncan_sex, pop_2_Duncan_age, pop_2_Duncan_educ)

pop_2_all_Duncan_round <- pop_2_all_Duncan %>% mutate(across(starts_with("Duncan"), round, 3))

Build table

pop_all_Duncan_round <- bind_cols(pop_0_all_Duncan_round, pop_1_all_Duncan_round, pop_2_all_Duncan_round)
pop_all_Duncan_round 

Moran’s index of spatial autocorrelation (Figure 1 and Appendix 3)

In night cells

dfpop_0_prop <- dfpop_0/dfpop_0$pop
pop_0_geom <- pop_0[, -1:-28] 
pop_0_prop <- bind_cols(pop_0_geom, dfpop_0_prop)

nb0 <- poly2nb(pop_0_prop, queen=TRUE)
lw0 <- nb2listw(nb0, style="W", zero.policy=TRUE)

list <- colnames(dfpop_0_prop[,2:27])
pop_0_moran <- data.frame(matrix(ncol = 2))
colnames(pop_0_moran) <- c("pop_night_IMoran" , "pop_night_Moranpvalue")

for (i in list)
{
  pop_0_i <- pop_0_prop[,i]
  colnames(pop_0_i)[1] ="ok"
  pop_0_moran[i, 1] <- moran.test(pop_0_i$ok, lw0, zero.policy=TRUE, alternative="greater") [3]
  pop_0_moran[i, 2] <- moran.test(pop_0_i$ok, lw0, zero.policy=TRUE, alternative="greater") [2]
}

pop_0_moran <- pop_0_moran[-1,] # Delete empty first line

pop_0_moran_round <- pop_0_moran %>% mutate(across(starts_with("pop"), round, 3))
pop_0_moran_round

In day cells

dfpop_1_prop <- dfpop_1/dfpop_1$pop
pop_1_geom <- pop_1[, -1:-28] 
pop_1_prop <- bind_cols(pop_1_geom, dfpop_1_prop)

nb1 <- poly2nb(pop_1_prop, queen=TRUE)
lw1 <- nb2listw(nb1, style="W", zero.policy=TRUE)

list <- colnames(dfpop_1_prop[,2:27])

pop_1_moran <- data.frame(matrix(ncol = 2))
colnames(pop_1_moran) <- c( "pop_day_IMoran" , "pop_day_Moranpvalue")

for (i in list)
{
  pop_1_i <- pop_1_prop[,i]
  colnames(pop_1_i)[1] ="ok"
  pop_1_moran[i, 1] <- moran.test(pop_1_i$ok, lw1, zero.policy=TRUE, alternative="greater") [3]
  pop_1_moran[i, 2] <- moran.test(pop_1_i$ok, lw1, zero.policy=TRUE, alternative="greater") [2]
}

pop_1_moran <- pop_1_moran[-1,] # Delete empty first line

pop_1_moran_round <- pop_1_moran %>% mutate(across(starts_with("pop"), round, 3))
pop_1_moran_round

In evening cells

dfpop_2_prop <- dfpop_2/dfpop_2$pop
pop_2_geom <- pop_2[, -1:-28] 
pop_2_prop <- bind_cols(pop_2_geom, dfpop_2_prop)

nb2 <- poly2nb(pop_2_prop, queen=TRUE)
lw2 <- nb2listw(nb2, style="W", zero.policy=TRUE)

list <- colnames(dfpop_2_prop[,2:27])

pop_2_moran <- data.frame(matrix(ncol = 2))
colnames(pop_2_moran) <- c( "pop_evening_IMoran" , "pop_evening_Moranpvalue")

for (i in list)
{
  pop_2_i <- pop_2_prop[,i]
  colnames(pop_2_i)[1] ="ok"
  pop_2_moran[i, 1] <- moran.test(pop_2_i$ok, lw2, zero.policy=TRUE, alternative="greater") [3]
  pop_2_moran[i, 2] <- moran.test(pop_2_i$ok, lw2, zero.policy=TRUE, alternative="greater") [2]
}

pop_2_moran <- pop_2_moran[-1,] # Delete empty first line

pop_2_moran_round <- pop_2_moran %>% mutate(across(starts_with("pop"), round, 3))
pop_2_moran_round

Combine Moran tables pop_0, pop_1, pop_2

pop_all_moran <- bind_cols(pop_0_moran_round, pop_1_moran_round, pop_2_moran_round)

Analysis of empirical data about dietary habits (section 2.2)

Summary by sociodemographic groups (cf. table 1)

data <- read.csv("data/sociodemo_distribution_healthy.csv", stringsAsFactors = F)
head(data)
##   Sex Age Edu   n_idf N_2002 N_2008 conso_5_2002 conso_5_2008
## 1   1   1   1  493871     52    148   0.01923077   0.02702703
## 2   1   1   2  496964     97    223   0.03092784   0.04932735
## 3   1   1   3  199476     84    137   0.03571429   0.07299270
## 4   1   2   1 1077670    241    355   0.04979253   0.08169014
## 5   1   2   2  662407    122    225   0.07377049   0.11111111
## 6   1   2   3  658465    118    180   0.05932203   0.16666667

n_idf = number of individuals of a particular group in 2012 in the Paris Region (Ile de France). N_2002 = number of individuals of a particular group in the Health and Nutrition Barometer survey in 2002. N_2008 = number of individuals of a particular group in the Health and Nutrition Barometer in 2008. conso_5_2002 = proportion of individuals eating 5+ portions of fruit and vegetables a day in the Health and Nutrition Barometer in 2002. conso_5_2008 = proportion of individuals eating 5+ portions of fruit and vegetables a day in the Health and Nutrition Barometer in 2008.

Computation of the total proportion of healthy individuals (i.e. eating 5+ portions of fruit and vegetables a day)

data$conso_5_2002 <- as.numeric(data$conso_5_2002)
data$conso_5_2008 <- as.numeric(data$conso_5_2008)
n_tot_idf <- sum(data$n_idf)
data$share_idf <- data$n_idf / n_tot_idf
data$h_idf_2002 <- round(data$n_idf * data$conso_5_2002,0)
data$h_idf_2008 <- round(data$n_idf * data$conso_5_2008,0)
data$h_baro_2002 <- round(data$N_2002 * data$conso_5_2002,0)
data$h_baro_2008 <- round(data$N_2008 * data$conso_5_2008,0)


total_prop <- function(data, year, sample){
  if(sample == "idf"){
    tp <- 
      sum(data[,paste("h", sample, year, sep="_")]) / 
      sum(data$n_idf)
  } else {
    tp <- 
      sum(data[,paste("h", sample, year, sep="_")]) / 
      sum(data[,paste("N", year, sep="_")])
  }
    return(tp)
  }

Results with the Paris region as reference population

total_prop(data=data, year = 2002, sample = "idf")
## [1] 0.09568032
total_prop(data=data, year = 2008, sample = "idf")
## [1] 0.1205141

Results with the diet survey as reference population

total_prop(data=data, year = 2002, sample = "baro")
## [1] 0.1063931
total_prop(data=data, year = 2008, sample = "baro")
## [1] 0.1201574

Computation of EducIneqIndex

n_tot_baro_2002 <- sum(data$N_2002)
n_tot_baro_2008 <- sum(data$N_2008)


EII <- function(data, year, sample){
  if(sample == "idf"){
    eii <- (
      data[data$Sex == 1 & data$Age == 1 & data$Edu == 3,
         paste0("conso_5_", year)] / 
       data[data$Sex == 1 & data$Age == 1 & data$Edu == 1,
            paste0("conso_5_", year)]
     ) * (
              sum(data[data$Sex == 1 & data$Age == 1, "n_idf"] / n_tot_idf) 
              ) + (
                data[data$Sex == 1 & data$Age == 2 & data$Edu == 3,
                     paste0("conso_5_", year)] / 
                  data[data$Sex == 1 & data$Age == 2 & data$Edu == 1,
                       paste0("conso_5_", year)]
              ) * (
                sum(data[data$Sex == 1 & data$Age == 2, "n_idf"] / n_tot_idf) 
              ) + (
                data[data$Sex == 1 & data$Age == 3 & data$Edu == 3,
                     paste0("conso_5_", year)] / 
                  data[data$Sex == 1 & data$Age == 3 & data$Edu == 1,
                       paste0("conso_5_", year)]
              ) * (
                sum(data[data$Sex == 1 & data$Age == 3, "n_idf"] / n_tot_idf) 
              ) + (
                data[data$Sex == 2 & data$Age == 1 & data$Edu == 3,
                     paste0("conso_5_", year)] / 
                  data[data$Sex == 2 & data$Age == 1 & data$Edu == 1,
                       paste0("conso_5_", year)]
              ) * (
                sum(data[data$Sex == 2 & data$Age == 1, "n_idf"] / n_tot_idf) 
              ) + (
                data[data$Sex == 2 & data$Age == 2 & data$Edu == 3,
                     paste0("conso_5_", year)] / 
                  data[data$Sex == 2 & data$Age == 2 & data$Edu == 1,
                       paste0("conso_5_", year)]
              ) * (
                sum(data[data$Sex == 2 & data$Age == 2, "n_idf"] / n_tot_idf) 
              ) + (
                data[data$Sex == 2 & data$Age == 3 & data$Edu == 3,
                     paste0("conso_5_", year)] / 
                  data[data$Sex == 2 & data$Age == 3 & data$Edu == 1,
                       paste0("conso_5_", year)]
              ) * (
                sum(data[data$Sex == 2 & data$Age == 3, "n_idf"] / n_tot_idf) 
              )
  } else {
    
    if(year == 2002){ n_tot_baro <- n_tot_baro_2002 }
    if(year == 2008){ n_tot_baro <- n_tot_baro_2008 }
    
    eii <- (
      data[data$Sex == 1 & data$Age == 1 & data$Edu == 3,
           paste0("conso_5_", year)] / 
        data[data$Sex == 1 & data$Age == 1 & data$Edu == 1,
             paste0("conso_5_", year)]
    ) * (
      sum(data[data$Sex == 1 & data$Age == 1, paste0("N_", year)] / n_tot_baro) 
    ) + (
      data[data$Sex == 1 & data$Age == 2 & data$Edu == 3,
           paste0("conso_5_", year)] / 
        data[data$Sex == 1 & data$Age == 2 & data$Edu == 1,
             paste0("conso_5_", year)]
    ) * (
      sum(data[data$Sex == 1 & data$Age == 2, paste0("N_", year)] / n_tot_baro) 
    ) + (
      data[data$Sex == 1 & data$Age == 3 & data$Edu == 3,
           paste0("conso_5_", year)] / 
        data[data$Sex == 1 & data$Age == 3 & data$Edu == 1,
             paste0("conso_5_", year)]
    ) * (
      sum(data[data$Sex == 1 & data$Age == 3, paste0("N_", year)] / n_tot_baro) 
    ) + (
      data[data$Sex == 2 & data$Age == 1 & data$Edu == 3,
           paste0("conso_5_", year)] / 
        data[data$Sex == 2 & data$Age == 1 & data$Edu == 1,
             paste0("conso_5_", year)]
    ) * (
      sum(data[data$Sex == 2 & data$Age == 1, paste0("N_", year)] / n_tot_baro) 
    ) + (
      data[data$Sex == 2 & data$Age == 2 & data$Edu == 3,
           paste0("conso_5_", year)] / 
        data[data$Sex == 2 & data$Age == 2 & data$Edu == 1,
             paste0("conso_5_", year)]
    ) * (
      sum(data[data$Sex == 2 & data$Age == 2, paste0("N_", year)] / n_tot_baro) 
    ) + (
      data[data$Sex == 2 & data$Age == 3 & data$Edu == 3,
           paste0("conso_5_", year)] / 
        data[data$Sex == 2 & data$Age == 3 & data$Edu == 1,
             paste0("conso_5_", year)]
    ) * (
      sum(data[data$Sex == 2 & data$Age == 3, paste0("N_", year)] / n_tot_baro) 
    )
      }
  return(eii)
}

Inequality results with the Paris region as reference population

EII(data=data, year = 2002, sample = "idf")
## [1] 1.407512
EII(data=data, year = 2008, sample = "idf")
## [1] 1.89476

Inequality results with the dietary survey as reference population

EII(data=data, year = 2002, sample = "baro")
## [1] 1.382566
EII(data=data, year = 2008, sample = "baro")
## [1] 1.905816

Evaluation of 5aday ABM model (section 3.4)

Importing calibration results

pareto <- read.csv2("data/nsga2_clean.csv")
head(pareto)
##   evolution.generation evolution.evaluated evolution.samples maxProbaToSwitch
## 1               376800              376800                 1            0.751
## 2               376800              376800                 6            0.750
## 3               376800              376800                 2            0.751
## 4               376800              376800                 9            0.749
## 5               376800              376800                 1            0.749
## 6               376800              376800                 1            0.748
##   constraintsStrength inertiaCoefficient healthyDietReward
## 1               0.100                  1                 0
## 2               0.100                  1                 0
## 3               0.098                  1                 0
## 4               0.097                  1                 0
## 5               0.097                  1                 0
## 6               0.097                  1                 0
##   objective.deltaHealth objective.deltaSocialInequality X.Healthy_vs._2002
## 1           331,736,427                           0.138               down
## 2           333,771,427                           0.149               down
## 3           334,153,927                           0.140               down
## 4           334,658,427                           0.128               down
## 5           335,441,427                           0.123               down
## 6           335,744,427                           0.123               down
##   X.Healthy_vs._2008 socialIneq_vs._2002         typo        typo_name
## 1               down                  up down_down_up decrease%Healthy
## 2               down                  up down_down_up decrease%Healthy
## 3               down                  up down_down_up decrease%Healthy
## 4               down                  up down_down_up decrease%Healthy
## 5               down                  up down_down_up decrease%Healthy
## 6               down                  up down_down_up decrease%Healthy
##                                                                                                                                                             socialInequality
## 1                                                                                                                                                       [1.7569191058549687]
## 2                                                          [1.7525945924779105,1.7540192352458608,1.7427673698247963,1.749312240252607,1.7412857735082057,1.742133022717344]
## 3                                                                                                                                     [1.7625539632685843,1.745986394640198]
## 4 [1.7662218959212503,1.7696432627743401,1.7665457668386697,1.756434185222282,1.7667110959021646,1.7696451049204427,1.7694793140516023,1.763192441014371,1.7589183006338713]
## 5                                                                                                                                                       [1.7713301959795948]
## 6                                                                                                                                                       [1.7720880068092864]
##                                                                                                                                                               deltaSocialInequality
## 1                                                                                                                                                             [0.13782823564503133]
## 2                                                         [0.14215274902208952,0.14072810625413923,0.15197997167520372,0.14543510124739312,0.15346156799179433,0.15261431878265608]
## 3                                                                                                                                         [0.13219337823141575,0.14876094685980212]
## 4 [0.12852544557874968,0.12510407872565987,0.1282015746613303,0.138313156277718,0.12803624559783544,0.12510223657955732,0.1252680274483977,0.13155490048562912,0.13582904086612868]
## 5                                                                                                                                                             [0.12341714552040517]
## 6                                                                                                                                                             [0.12265933469071366]
##                                                                                                            opinionMedian
## 1                                                                                                  [0.49988672137260437]
## 2 [0.4999333918094635,0.4999331533908844,0.49993330240249634,0.49993327260017395,0.4999333322048187,0.49993306398391724]
## 3                                                                               [0.4997502863407135,0.49975064396858215]
## 4                                                                                  [0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5]
## 5                                                                                                                  [0.5]
## 6                                                                                                    [0.499999463558197]
##                                                    numberOfHealthy
## 1                                                         [834951]
## 2                      [834948,832940,835222,834168,833731,834948]
## 3                                                  [827118,826925]
## 4 [818522,817720,816927,815986,817833,817780,815267,817438,818191]
## 5                                                         [817706]
## 6                                                         [811529]
##                                                                                                                                                                    deltaHealth
## 1                                                                                                                                                         [331736.42711645283]
## 2                                                           [332013.5408605285,334701.42711645283,332841.42711645283,332113.42711645283,335052.42711645283,334725.42711645283]
## 3                                                                                                                                      [333369.42711645283,334938.42711645283]
## 4 [333959.42711645283,333931.42711645283,335330.42711645283,334215.42711645283,336478.42711645283,336057.42711645283,334658.42711645283,335249.42711645283,334452.42711645283]
## 5                                                                                                                                                         [335441.42711645283]
## 6                                                                                                                                                         [335744.42711645283]

Drawing the Pareto front to select a parameter set (figure 4)

pareto$reliability.sample <- ifelse(pareto$evolution.samples>5, 2,1)
pareto$objective.deltaHealth <- as.numeric(gsub(",","",pareto$objective.deltaHealth))

front <- ggplot(data=pareto, aes(x = objective.deltaHealth, y=objective.deltaSocialInequality))
front + 
  geom_point(aes(group = typo_name, colour = typo_name, 
                       size=reliability.sample)) +
  scale_size_continuous(range = c(0.5,1.5)) + 
  scale_colour_manual(values = alpha(
    c("#faa911", #orange
      "#07b55e", #green
      "#1088b0", #blue
      "#aae860" #light green
      ), .5)) +
  geom_vline(xintercept = 290000000) + 
  geom_hline(yintercept = 0.36) +
  theme_bw() 

Drawing more figures about the Pareto front (Appendix 6)

pareto$id <- rownames(pareto)
colnames(pareto)
##  [1] "evolution.generation"            "evolution.evaluated"            
##  [3] "evolution.samples"               "maxProbaToSwitch"               
##  [5] "constraintsStrength"             "inertiaCoefficient"             
##  [7] "healthyDietReward"               "objective.deltaHealth"          
##  [9] "objective.deltaSocialInequality" "X.Healthy_vs._2002"             
## [11] "X.Healthy_vs._2008"              "socialIneq_vs._2002"            
## [13] "typo"                            "typo_name"                      
## [15] "socialInequality"                "deltaSocialInequality"          
## [17] "opinionMedian"                   "numberOfHealthy"                
## [19] "deltaHealth"                     "reliability.sample"             
## [21] "id"
pareto_long <- melt(pareto[,c("maxProbaToSwitch","constraintsStrength" ,
                              "inertiaCoefficient","healthyDietReward",
                              "objective.deltaSocialInequality", "objective.deltaHealth",
                              "typo_name", "reliability.sample")], 
                    id.vars=c("objective.deltaSocialInequality", "objective.deltaHealth", 
                              "typo_name", "reliability.sample"))

ggplot(pareto_long, aes(objective.deltaHealth, value)) + 
  geom_point(aes(group = typo_name, colour = typo_name, 
                 size=reliability.sample)) +
  scale_size_continuous(range = c(0.5,1.5)) + scale_colour_manual(values = alpha(
    c("#faa911", #orange
      "#07b55e", #green
      "#1088b0", #blue
      "#aae860" #light green
    ), .5)) +
  theme_bw() +  theme(legend.position="none") +
  facet_wrap(~variable, scales = "free_y", ncol = 1)

ggplot(pareto_long, aes(value, objective.deltaSocialInequality)) + 
  geom_point(aes(group = typo_name, colour = typo_name, 
                 size=reliability.sample)) +
  scale_size_continuous(range = c(0.5,1.5)) + scale_colour_manual(values = alpha(
    c("#faa911", #orange
      "#07b55e", #green
      "#1088b0", #blue
      "#aae860" #light green
    ), .5)) +
  theme_bw()  +  theme(legend.position="none")+
  facet_wrap(~variable, scales = "free_x", nrow = 1)

Results (section 4)

Social inequalities in healthy behaviour for the 5 scenarios (section 4.1)

Importing simulation results

sim <- read.csv2("data/results-25042022-10000replications.csv", sep=",", dec=".")
sim$Scenario <- ifelse(sim$numOfScenario == 1, "1A: Random residence / no move",
                       ifelse(sim$numOfScenario == 2, "1B: Random residence / random moves",
                              ifelse(sim$numOfScenario == 3, "2A: Observed residence / no move",
                                     ifelse(sim$numOfScenario == 4, "2B: Observed residence / random moves",
                                            "2C: Observed residence / observed moves"))))
head(sim)
##   constraintsStrength erreygersE healthyDietReward inertiaCoefficient
## 1            0.128944 0.03559005          0.153103          0.8195314
## 2            0.128944 0.03589710          0.153103          0.8195314
## 3            0.128944 0.03558682          0.153103          0.8195314
## 4            0.128944 0.03597858          0.153103          0.8195314
## 5            0.128944 0.03580784          0.153103          0.8195314
## 6            0.128944 0.03594730          0.153103          0.8195314
##   maxProbaToSwitch numOfScenario numberOfHealthy          seed socialInequality
## 1        0.8904704             1         1187441  5.241114e+18         1.384503
## 2        0.8904704             1         1188705  8.182251e+18         1.389434
## 3        0.8904704             1         1188073 -5.732413e+18         1.383992
## 4        0.8904704             1         1190153 -9.124344e+18         1.388518
## 5        0.8904704             1         1189071 -3.952058e+18         1.385801
## 6        0.8904704             1         1191965 -7.961473e+18         1.385447
##   openmole.experiment                       Scenario
## 1               20103 1A: Random residence / no move
## 2               20099 1A: Random residence / no move
## 3               20107 1A: Random residence / no move
## 4               20105 1A: Random residence / no move
## 5               20101 1A: Random residence / no move
## 6               20285 1A: Random residence / no move

Distribution of simulated values of EducIneqIndex (EII) for the five space-time scenarios (10 000 replications per scenario) (figure 5)

mean_1A <- mean(sim[sim$numOfScenario == 1,"socialInequality"])
mean_1B <- mean(sim[sim$numOfScenario == 2,"socialInequality"])
mean_2A <- mean(sim[sim$numOfScenario == 3,"socialInequality"])
mean_2B <- mean(sim[sim$numOfScenario == 4,"socialInequality"])
mean_2C <- mean(sim[sim$numOfScenario == 5,"socialInequality"])

sim %>%
  ggplot( aes(x=socialInequality, group=Scenario, fill=Scenario)) +
  geom_vline(xintercept=1.405, size=1, color="grey30", linetype = "dashed") +
  geom_density(alpha=0.6) +
  scale_fill_manual(values=c("firebrick2", "mediumorchid4", "dodgerblue2", "darkorange2", "springgreen3"))+
  theme_bw() +
  geom_vline(xintercept=mean_1A, size=0.4, color="firebrick2") +
  geom_vline(xintercept=mean_1B, size=0.4, color="mediumorchid4") +
  geom_vline(xintercept=mean_2A, size=0.4, color="dodgerblue2") +
  geom_vline(xintercept=mean_2B, size=0.4, color="darkorange2") +
  geom_vline(xintercept=mean_2C, size=0.4, color="springgreen3") +
  annotate("text", x=mean_1A - 0.004, y= -10, size=2.5, label=round(mean_1A,3), color="firebrick2") +
  annotate("text", x=mean_1B - 0.004 , y= -10, size=2.5, label=round(mean_1B,3), color="mediumorchid4") +
  annotate("text", x=mean_2A - 0.004 , y= -10, size=2.5, label=round(mean_2A,3), color="dodgerblue2") +
  annotate("text", x=mean_2B + 0.004 , y= -10, size=2.5, label=round(mean_2B,3), color="darkorange2") +
  annotate("text", x=mean_2C + 0.004, y= -10, size=2.5, label=round(mean_2C,3), color="springgreen3") +
  annotate("text", x=1.405 + 0.017, y= 115, size=3, label="Value at initialisation", color="grey30") +
  xlab("Social Inequality Index") +
  ylab("Number of simulations") + 
  guides(fill=guide_legend(title="Type of scenario")) #+ 

Overall proportions of agents with healthy behaviour for the 5 scenarios (section 4.2)

!! TO DO JULIE (Table 5)

Spatial distribution of agents with healthy behaviour for the 5 scenarios - in the night cells (section 4.3)

s0r <- st_read( "data/results_HigherProp_Scenario1_RandomPop_NoMove_42.gpkg",
              layer="Initial_0", quiet = TRUE)
s0o <- st_read( "data/results_HigherProp_Scenario5_ObservedPop_ObservedMove_42.gpkg",
               layer="Initial_0", quiet = TRUE)

s1A <- st_read( "data/results_HigherProp_Scenario1_RandomPop_NoMove_42.gpkg",
               layer="Final", quiet = TRUE)
s1B <- st_read( "data/results_HigherProp_Scenario2_RandomPop_RandomMove_42.gpkg",
               layer="Final", quiet = TRUE)

s2A <- st_read( "data/results_HigherProp_Scenario3_ObservedPop_NoMove_42.gpkg",
               layer="Final", quiet = TRUE)
s2B <- st_read( "data/results_HigherProp_Scenario4_ObservedPop_RandomMove_42.gpkg",
               layer="Final", quiet = TRUE)
s2C <- st_read( "data/results_HigherProp_Scenario5_ObservedPop_ObservedMove_42.gpkg",
               layer="Final", quiet = TRUE)

bks <- c(0, 0.1, 0.125, 0.15, 0.175, Inf)

Maps (figure 6)

bks <- c(0, 0.1, 0.125, 0.15, 0.175, Inf)

Initial Random

tm_shape(s0r) + tm_fill(col="X.propHealthy", style="fixed", 
                       breaks = bks, palette="Greens") +
  tm_legend(outside=TRUE)

Initial Observed

tm_shape(s0o) + tm_fill(col="X.propHealthy", style="fixed", 
                       breaks = bks, palette="Greens") +
  tm_legend(outside=TRUE)

Scenario 1A

tm_shape(s1A) + tm_fill(col="X.propHealthy", style="fixed", 
                       breaks = bks, palette="Greens") +
  tm_legend(outside=TRUE)

Scenario 1B

tm_shape(s1B) + tm_fill(col="X.propHealthy", style="fixed", 
                        breaks = bks, palette="Greens") +
  tm_legend(outside=TRUE)

Scenario 2A

tm_shape(s2A) + tm_fill(col="X.propHealthy", style="fixed", 
                        breaks = bks, palette="Greens") +
  tm_legend(outside=TRUE)

Scenario 2B

tm_shape(s2B) + tm_fill(col="X.propHealthy", style="fixed", 
                        breaks = bks, palette="Greens") +
  tm_legend(outside=TRUE)

Scenario 2C

tm_shape(s2C) + tm_fill(col="X.propHealthy", style="fixed", 
                        breaks = bks, palette="Greens") +
  tm_legend(outside=TRUE)

Duncan Dissimilarity Index (figure 6)

duncan <- data.frame(matrix(ncol = 2, nrow=7))
colnames(duncan) <- c("type", "Duncan" )

dfs0r<- s0r %>% st_drop_geometry()
duncan[1,1] <- "s0r_night"
duncan[1,2]<- ISDuncan(dfs0r[,3:4]) [1]

dfs0o<- s0o %>% st_drop_geometry()
duncan[2,1] <- "s0o_night"
duncan[2,2]<- ISDuncan(dfs0o[,3:4]) [1]


dfs1A<- s1A %>% st_drop_geometry()
duncan[3,1] <- "s1A_final_night"
duncan[3,2]<- ISDuncan(dfs1A[,3:4]) [1]

dfs1B<- s1B %>% st_drop_geometry()
duncan[4,1] <- "s1B_final_night"
duncan[4,2]<- ISDuncan(dfs1B[,3:4]) [1]


dfs2A<- s2A %>% st_drop_geometry()
duncan[5,1] <- "s2A_final_night"
duncan[5,2]<- ISDuncan(dfs2A[,3:4]) [1]

dfs2B<- s2B %>% st_drop_geometry()
duncan[6,1] <- "s2B_final_night"
duncan[6,2]<- ISDuncan(dfs2B[,3:4]) [1]

dfs2C<- s2C %>% st_drop_geometry()
duncan[7,1] <- "s2C_final_night"
duncan[7,2]<- ISDuncan(dfs2C[,3:4]) [1]

duncan_round <- duncan %>% mutate(across(starts_with("Duncan"), round, 3))

duncan_round

Moran’s index of spatial autocorrelation (figure 6)

moran <- data.frame(matrix(ncol = 3, nrow=7))
colnames(moran) <- c("type", "IMoran" , "pvalueMoran")

nbs0r <- poly2nb(s0r, queen=TRUE)
lws0r  <- nb2listw(nbs0r, style="W", zero.policy=TRUE)
moran[1,1] <- "s0r_night"
moran[1,2] <- moran.test(s0r$X.propHealthy, lws0r, zero.policy=TRUE, alternative="greater") [3]
moran[1,3] <- moran.test(s0r$X.propHealthy, lws0r, zero.policy=TRUE, alternative="greater") [2]

nbs0o <- poly2nb(s0o, queen=TRUE)
lws0o  <- nb2listw(nbs0o, style="W", zero.policy=TRUE)
moran[4,1] <- "s0o_night"
moran[4,2] <- moran.test(s0o$X.propHealthy, lws0o, zero.policy=TRUE, alternative="greater") [3]
moran[4,3] <- moran.test(s0o$X.propHealthy, lws0o, zero.policy=TRUE, alternative="greater") [2]

nbs1A <- poly2nb(s1A, queen=TRUE)
lws1A  <- nb2listw(nbs1A, style="W", zero.policy=TRUE)
moran[2,1] <- "s1A_final_night"
moran[2,2] <- moran.test(s1A$X.propHealthy, lws1A, zero.policy=TRUE, alternative="greater") [3]
moran[2,3] <- moran.test(s1A$X.propHealthy, lws1A, zero.policy=TRUE, alternative="greater") [2]

nbs1B <- poly2nb(s1B, queen=TRUE)
lws1B  <- nb2listw(nbs1B, style="W", zero.policy=TRUE)
moran[3,1] <- "s1B_final_night"
moran[3,2] <- moran.test(s1B$X.propHealthy, lws1B, zero.policy=TRUE, alternative="greater") [3]
moran[3,3] <- moran.test(s1B$X.propHealthy, lws1B, zero.policy=TRUE, alternative="greater") [2]


nbs2A <- poly2nb(s2A, queen=TRUE)
lws2A  <- nb2listw(nbs2A, style="W", zero.policy=TRUE)
moran[5,1] <- "s2A_final_night"
moran[5,2] <- moran.test(s2A$X.propHealthy, lws2A, zero.policy=TRUE, alternative="greater") [3]
moran[5,3] <- moran.test(s2A$X.propHealthy, lws2A, zero.policy=TRUE, alternative="greater") [2]

nbs2B <- poly2nb(s2B, queen=TRUE)
lws2B  <- nb2listw(nbs2B, style="W", zero.policy=TRUE)
moran[6,1] <- "s2B_final_night"
moran[6,2] <- moran.test(s2B$X.propHealthy, lws2B, zero.policy=TRUE, alternative="greater") [3]
moran[6,3] <- moran.test(s2B$X.propHealthy, lws2B, zero.policy=TRUE, alternative="greater") [2]

nbs2C <- poly2nb(s2C, queen=TRUE)
lws2C  <- nb2listw(nbs2C, style="W", zero.policy=TRUE)
moran[7,1] <- "s2C_final_night"
moran[7,2] <- moran.test(s2C$X.propHealthy, lws2C, zero.policy=TRUE, alternative="greater") [3]
moran[7,3] <- moran.test(s2C$X.propHealthy, lws2C, zero.policy=TRUE, alternative="greater") [2]

moran_round <- moran %>% mutate(across(2:3, round, 3))

moran_round 
##              type IMoran pvalueMoran
## 1       s0r_night  0.010       0.061
## 2 s1A_final_night  0.010       0.057
## 3 s1B_final_night  0.016       0.008
## 4       s0o_night  0.006       0.175
## 5 s2A_final_night  0.033       0.000
## 6 s2B_final_night  0.024       0.000
## 7 s2C_final_night  0.009       0.091

Alternative measure of social inequality in health behaviour (Appendix 4)

Computation of rank-ordered inequality measure

data_edu <- aggregate(data[,c("n_idf", "N_2002", "N_2008",
                              "h_idf_2002", "h_idf_2008", 
                              "h_baro_2002", "h_baro_2008")], 
                      by = list(data$Edu), FUN = sum)
colnames(data_edu)[1] <- "edu"
data_edu
##   edu   n_idf N_2002 N_2008 h_idf_2002 h_idf_2008 h_baro_2002 h_baro_2008
## 1   1 3800999   1030   1372     336079     359509         122         155
## 2   2 2946028    595   1135     279323     361858          57         132
## 3   3 2021871    471    797     223609     335409          44         110
guido <- function(data_edu, year, sample){
  
  if(sample == "idf"){
    n <- sum(data_edu[,"n_idf"])
    n1 <- data_edu[data_edu$edu == 1,"n_idf"]
    n2 <- data_edu[data_edu$edu == 2,"n_idf"]
    n3 <- data_edu[data_edu$edu == 3,"n_idf"]
  } else {
    n <- sum(data_edu[,paste("N", year, sep="_")])
    n1 <- data_edu[data_edu$edu == 1,paste("N", year, sep="_")]
    n2 <- data_edu[data_edu$edu == 2,paste("N", year, sep="_")]
    n3 <- data_edu[data_edu$edu == 3,paste("N", year, sep="_")]
  }
  
  h1 <- data_edu[data_edu$edu == 1,paste("h", sample, year, sep="_")]
  h2 <- data_edu[data_edu$edu == 2,paste("h", sample, year, sep="_")]
  h3 <- data_edu[data_edu$edu == 3,paste("h", sample, year, sep="_")]
  
  t1 <- (1 - n2 - n3) / 2
  t2 <- (n1 - n3 + 1) / 2
  t3 <- (n1 + n2 + 1) / 2
 
  guido <- ( 8 / n^2 ) * (h1 * t1 + h2 * t2 + h3 * t3)
  return(guido)
}

guido(data_edu = data_edu, year = 2002, sample = "idf")
## [1] 0.01748087
guido(data_edu = data_edu, year = 2008, sample = "idf")
## [1] 0.05830406
guido(data_edu = data_edu, year = 2002, sample = "baro")
## [1] -0.02409715
guido(data_edu = data_edu, year = 2008, sample = "baro")
## [1] 0.01927629

Density distribution of rank-order inequality index per scenario

mean_1A <- mean(sim[sim$numOfScenario == 1,"erreygersE"])
mean_1B <- mean(sim[sim$numOfScenario == 2,"erreygersE"])
mean_2A <- mean(sim[sim$numOfScenario == 3,"erreygersE"])
mean_2B <- mean(sim[sim$numOfScenario == 4,"erreygersE"])
mean_2C <- mean(sim[sim$numOfScenario == 5,"erreygersE"])

sim %>%
  ggplot( aes(x=erreygersE, group=Scenario, fill=Scenario)) +
  geom_vline(xintercept=0.0175, size=1, color="grey30", linetype = "dashed") +
  geom_density(alpha=0.6) +
  scale_fill_manual(values=c("firebrick2", "mediumorchid4", "dodgerblue2", "darkorange2", "springgreen3"))+
  theme_bw() +
  geom_vline(xintercept=mean_1A, size=0.4, color="firebrick2") +
  geom_vline(xintercept=mean_1B, size=0.4, color="mediumorchid4") +
  geom_vline(xintercept=mean_2A, size=0.4, color="dodgerblue2") +
  geom_vline(xintercept=mean_2B, size=0.4, color="darkorange2") +
  geom_vline(xintercept=mean_2C, size=0.4, color="springgreen3") +
  annotate("text", x=mean_1A - 0.004, y= -10, size=2.5, label=round(mean_1A,3), color="firebrick2") +
  annotate("text", x=mean_1B - 0.004 , y= -10, size=2.5, label=round(mean_1B,3), color="mediumorchid4") +
  annotate("text", x=mean_2A - 0.004 , y= -10, size=2.5, label=round(mean_2A,3), color="dodgerblue2") +
  annotate("text", x=mean_2B + 0.004 , y= -10, size=2.5, label=round(mean_2B,3), color="darkorange2") +
  annotate("text", x=mean_2C + 0.004, y= -10, size=2.5, label=round(mean_2C,3), color="springgreen3") +
  annotate("text", x=0.0175 + 0.0051, y= 115, size=3, label="Value at initialisation", color="grey30") +
  xlab("Rank-Order Inequality Index") +
  ylab("Number of simulations") + 
  guides(fill=guide_legend(title="Type of scenario")) #+ 

Alternative measures of spatial distribution of agents with healthy behaviour - in the day and evening cells (Appendix 8)

In day cells

s1A_day <- st_read( "data/results_HigherProp_Scenario1_RandomPop_NoMove_42.gpkg", quiet = TRUE,
                layer="Final_1")

s1B_day <- st_read( "data/results_HigherProp_Scenario2_RandomPop_RandomMove_42.gpkg", quiet = TRUE,
                layer="Final_1")

s2A_day <- st_read( "data/results_HigherProp_Scenario3_ObservedPop_NoMove_42.gpkg", quiet = TRUE,
                layer="Final_1")

s2B_day <- st_read( "data/results_HigherProp_Scenario4_ObservedPop_RandomMove_42.gpkg", quiet = TRUE,
                layer="Final_1")

s2C_day <- st_read( "data/results_HigherProp_Scenario5_ObservedPop_ObservedMove_42.gpkg", quiet = TRUE,
                layer="Final_1")

Duncan

duncan_day <- data.frame(matrix(ncol = 2, nrow=5))
colnames(duncan_day) <- c("type", "Duncan_dayCell" )

dfs1A_day<- s1A_day %>% st_drop_geometry()
duncan_day[1,1] <- "s1A_final_day"
duncan_day[1,2]<- ISDuncan(dfs1A_day[,3:4]) [1]

dfs1B_day<- s1B_day %>% st_drop_geometry()
duncan_day[2,1] <- "s1B_final_day"
duncan_day[2,2]<- ISDuncan(dfs1B_day[,3:4]) [1]

dfs2A_day<- s2A_day %>% st_drop_geometry()
duncan_day[3,1] <- "s2A_final_day"
duncan_day[3,2]<- ISDuncan(dfs2A_day[,3:4]) [1]

dfs2B_day<- s2B_day %>% st_drop_geometry()
duncan_day[4,1] <- "s2B_final_day"
duncan_day[4,2]<- ISDuncan(dfs2B_day[,3:4]) [1]

dfs2C_day<- s2C_day %>% st_drop_geometry()
duncan_day[5,1] <- "s2C_final_day"
duncan_day[5,2]<- ISDuncan(dfs2C_day[,3:4]) [1]

duncan_day <- duncan_day %>% mutate(across(starts_with("Duncan"), round, 3))

Moran

moran_day <- data.frame(matrix(ncol = 3, nrow=5))
colnames(moran_day) <- c("type", "IMoran_dayCell" , "pvalueMoran_dayCell")

nbs1A_day <- poly2nb(s1A_day, queen=TRUE)
lws1A_day  <- nb2listw(nbs1A_day, style="W", zero.policy=TRUE)
moran_day[1,1] <- "s1A_final_day"
moran_day[1,2] <- moran.test(s1A_day$X.propHealthy, lws1A_day, zero.policy=TRUE, alternative="greater") [3]
moran_day[1,3] <- moran.test(s1A_day$X.propHealthy, lws1A_day, zero.policy=TRUE, alternative="greater") [2]

nbs1B_day <- poly2nb(s1B_day, queen=TRUE)
lws1B_day  <- nb2listw(nbs1B_day, style="W", zero.policy=TRUE)
moran_day[2,1] <- "s1B_final_day"
moran_day[2,2] <- moran.test(s1B_day$X.propHealthy, lws1B_day, zero.policy=TRUE, alternative="greater") [3]
moran_day[2,3] <- moran.test(s1B_day$X.propHealthy, lws1B_day, zero.policy=TRUE, alternative="greater") [2]

nbs2A_day <- poly2nb(s2A_day, queen=TRUE)
lws2A_day  <- nb2listw(nbs2A_day, style="W", zero.policy=TRUE)
moran_day[3,1] <- "s2A_final_day"
moran_day[3,2] <- moran.test(s2A_day$X.propHealthy, lws2A_day, zero.policy=TRUE, alternative="greater") [3]
moran_day[3,3] <- moran.test(s2A_day$X.propHealthy, lws2A_day, zero.policy=TRUE, alternative="greater") [2]

nbs2B_day <- poly2nb(s2B_day, queen=TRUE)
lws2B_day  <- nb2listw(nbs2B_day, style="W", zero.policy=TRUE)
moran_day[4,1] <- "s2B_final_day"
moran_day[4,2] <- moran.test(s2B_day$X.propHealthy, lws2B_day, zero.policy=TRUE, alternative="greater") [3]
moran_day[4,3] <- moran.test(s2B_day$X.propHealthy, lws2B_day, zero.policy=TRUE, alternative="greater") [2]

nbs2C_day <- poly2nb(s2C_day, queen=TRUE)
lws2C_day  <- nb2listw(nbs2C_day, style="W", zero.policy=TRUE)
moran_day[5,1] <- "s2C_final_day"
moran_day[5,2] <- moran.test(s2C_day$X.propHealthy, lws2C_day, zero.policy=TRUE, alternative="greater") [3]
moran_day[5,3] <- moran.test(s2C_day$X.propHealthy, lws2C_day, zero.policy=TRUE, alternative="greater") [2]

moran_round_day <- moran_day %>% mutate(across(2:3, round, 3))

Combine Duncan and Moran files

Duncan_moran_dayCell <- left_join(duncan_day, moran_round_day, by = "type")
Duncan_moran_dayCell 

In evening cells

s1A_evening <- st_read( "data/results_HigherProp_Scenario1_RandomPop_NoMove_42.gpkg", quiet = TRUE,
                    layer="Final_2")

s1B_evening <- st_read( "data/results_HigherProp_Scenario2_RandomPop_RandomMove_42.gpkg", quiet = TRUE,
                    layer="Final_2")

s2A_evening <- st_read( "data/results_HigherProp_Scenario3_ObservedPop_NoMove_42.gpkg", quiet = TRUE,
                    layer="Final_2")

s2B_evening <- st_read( "data/results_HigherProp_Scenario4_ObservedPop_RandomMove_42.gpkg", quiet = TRUE,
                    layer="Final_2")

s2C_evening <- st_read( "data/results_HigherProp_Scenario5_ObservedPop_ObservedMove_42.gpkg", quiet = TRUE,
                    layer="Final_2")

Duncan

duncan_evening <- data.frame(matrix(ncol = 2, nrow=5))
colnames(duncan_evening) <- c("type", "Duncan_eveningCell" )


dfs1A_evening<- s1A_evening %>% st_drop_geometry()
duncan_evening[1,1] <- "s1A_final_evening"
duncan_evening[1,2]<- ISDuncan(dfs1A_evening[,3:4]) [1]

dfs1B_evening<- s1B_evening %>% st_drop_geometry()
duncan_evening[2,1] <- "s1B_final_evening"
duncan_evening[2,2]<- ISDuncan(dfs1B_evening[,3:4]) [1]

dfs2A_evening<- s2A_evening %>% st_drop_geometry()
duncan_evening[3,1] <- "s2A_final_evening"
duncan_evening[3,2]<- ISDuncan(dfs2A_evening[,3:4]) [1]

dfs2B_evening<- s2B_evening %>% st_drop_geometry()
duncan_evening[4,1] <- "s2B_final_evening"
duncan_evening[4,2]<- ISDuncan(dfs2B_evening[,3:4]) [1]

dfs2C_evening<- s2C_evening %>% st_drop_geometry()
duncan_evening[5,1] <- "s2C_final_evening"
duncan_evening[5,2]<- ISDuncan(dfs2C_evening[,3:4]) [1]

duncan_evening <- duncan_evening %>% mutate(across(starts_with("Duncan"), round, 3))

Moran

moran_evening <- data.frame(matrix(ncol = 3, nrow=5))
colnames(moran_evening) <- c("type", "IMoran_eveningCell" , "pvalueMoran_eveningCell")

nbs1A_evening <- poly2nb(s1A_evening, queen=TRUE)
lws1A_evening  <- nb2listw(nbs1A_evening, style="W", zero.policy=TRUE)
moran_evening[1,1] <- "s1A_final_evening"
moran_evening[1,2] <- moran.test(s1A_evening$X.propHealthy, lws1A_evening, zero.policy=TRUE, alternative="greater") [3]
moran_evening[1,3] <- moran.test(s1A_evening$X.propHealthy, lws1A_evening, zero.policy=TRUE, alternative="greater") [2]

nbs1B_evening <- poly2nb(s1B_evening, queen=TRUE)
lws1B_evening  <- nb2listw(nbs1B_evening, style="W", zero.policy=TRUE)
moran_evening[2,1] <- "s1B_final_evening"
moran_evening[2,2] <- moran.test(s1B_evening$X.propHealthy, lws1B_evening, zero.policy=TRUE, alternative="greater") [3]
moran_evening[2,3] <- moran.test(s1B_evening$X.propHealthy, lws1B_evening, zero.policy=TRUE, alternative="greater") [2]

nbs2A_evening <- poly2nb(s2A_evening, queen=TRUE)
lws2A_evening  <- nb2listw(nbs2A_evening, style="W", zero.policy=TRUE)
moran_evening[3,1] <- "s2A_final_evening"
moran_evening[3,2] <- moran.test(s2A_evening$X.propHealthy, lws2A_evening, zero.policy=TRUE, alternative="greater") [3]
moran_evening[3,3] <- moran.test(s2A_evening$X.propHealthy, lws2A_evening, zero.policy=TRUE, alternative="greater") [2]

nbs2B_evening <- poly2nb(s2B_evening, queen=TRUE)
lws2B_evening  <- nb2listw(nbs2B_evening, style="W", zero.policy=TRUE)
moran_evening[4,1] <- "s2B_final_evening"
moran_evening[4,2] <- moran.test(s2B_evening$X.propHealthy, lws2B_evening, zero.policy=TRUE, alternative="greater") [3]
moran_evening[4,3] <- moran.test(s2B_evening$X.propHealthy, lws2B_evening, zero.policy=TRUE, alternative="greater") [2]

nbs2C_evening <- poly2nb(s2C_evening, queen=TRUE)
lws2C_evening  <- nb2listw(nbs2C_evening, style="W", zero.policy=TRUE)
moran_evening[5,1] <- "s2C_final_evening"
moran_evening[5,2] <- moran.test(s2C_evening$X.propHealthy, lws2C_evening, zero.policy=TRUE, alternative="greater") [3]
moran_evening[5,3] <- moran.test(s2C_evening$X.propHealthy, lws2C_evening, zero.policy=TRUE, alternative="greater") [2]

moran_round_evening <- moran_evening %>% mutate(across(2:3, round, 3))

Combine Duncan and Moran files

Duncan_moran_eveningCell <- left_join(duncan_evening, moran_round_evening, by = "type")
Duncan_moran_eveningCell

Evolution over the six series of three time slices in simulated values - 5aDay model (Appendix 9)

library(cowplot)
library(tidyverse)
library(here)
 specialRead <- function(x) {
   read_csv(file=x, col_names = TRUE, col_types= "nnnnnnnnnn")
 }

concatFiles <- function (base, folder){
   fullPath = paste(base,folder,sep="")
   list.files(path = fullPath, full.names = TRUE, pattern = "(\\d_\\d_\\d)") %>% lapply(specialRead) %>% bind_rows
}

recode <- function(x) {
  step_to_day <- c (
  '00' = 'day 0 step 1',
  '01' = 'day 0 step 2',
  '02' = 'day 0 step 3',
  '10' = 'day 1 step 1',
  '11' = 'day 1 step 2',
  '12' = 'day 1 step 3',
  '20' = 'day 2 step 1',
  '21' = 'day 2 step 2',
  '22' = 'day 2 step 3',
  '30' = 'day 3 step 1',
  '31' = 'day 3 step 2',
  '32' = 'day 3 step 3',
  '40' = 'day 4 step 1',
  '41' = 'day 4 step 2',
  '42' = 'day 4 step 3',
  '50' = 'day 5 step 1',
  '51' = 'day 5 step 2',
  '52' = 'day 5 step 3'
  )
  return(step_to_day[x])
}

Simulated proportions of healthy agents, of social inequality (EII) and of rank-ordered inequalities (Et)

 hp_sc5_health <- read_csv("data/HigherProp/health.csv")
 hp_sc5_health$propHealthy <- hp_sc5_health$healthy / hp_sc5_health$effective

 hpLong_sc5_health <- hp_sc5_health %>%
   gather(key=facteur, value=valeur, propHealthy,avgOpinion,socialInequality,e ) %>%
   mutate(ds= recode(paste0(day, slice)))

 hpLong_sc5_health <-  hpLong_sc5_health %>% mutate (type="hp")

 both_sc5_health <- hpLong_sc5_health %>% mutate(type = as.factor(type))
 p_hp_propHealthy <- both_sc5_health %>% filter(facteur=="propHealthy")  %>% ggplot( aes(x=ds, y=valeur,colour=type) )+
   theme_bw()+
   theme(axis.text.x = element_text(angle = 45, vjust = 0.9, hjust=1.0)) +
   theme(axis.text.x.bottom = element_text(margin=margin(t=0, r=0,b=10,l=0)) ) +
   geom_point() +
    theme(legend.position = "none") +
   labs(y= "Value", x = "prop. of healthy")

 p_hp_ineq <- both_sc5_health %>% filter(facteur=="socialInequality")  %>% ggplot( aes(x=ds, y=valeur,colour=type)) +
   theme_bw()+
   theme(axis.text.x = element_text(angle = 45, vjust = 0.9, hjust=1.0)) +
   theme(axis.text.x.bottom = element_text(margin=margin(t=0, r=0,b=10,l=0)) ) +
   geom_point() +
   theme(legend.position = "none") +
   labs(y= "Value", x = "EducIneqIndex (EII)")

 p_hp_e <- both_sc5_health %>% filter(facteur=="e")  %>% ggplot( aes(x=ds, y=valeur,colour=type)) +
   theme_bw()+
   theme(axis.text.x = element_text(angle = 45, vjust = 0.9, hjust=1.0)) +
   theme(axis.text.x.bottom = element_text(margin=margin(t=0, r=0,b=10,l=0)) ) +
   geom_point() +
   theme(legend.position = "none") +
   labs(y= "Value", x = "Rank-ordered inequalities (Et)")

 p_hp_all <- plot_grid(p_hp_propHealthy, p_hp_ineq, p_hp_e, labels = "AUTO",nrow=1, align="h")
 p_hp_all

 ##ggsave("out/p_hp_all.pdf",p_hp_all, width=20)

Simulated proportion of healthy agents and simulated opinion across the 18 sociodemographic categories

sortie_sc5_higherprop <- concatFiles("data", "/HigherProp")
sortie_sc5_higherprop$propHealthy <- sortie_sc5_higherprop$healthy / sortie_sc5_higherprop$effective

sortieLong_sc5_hp <- gather(sortie_sc5_higherprop, facteur, valeur, c(propHealthy, avgOpinion) )
sortieLong_sc5_hp <- mutate(sortieLong_sc5_hp, "ds" = gsub(" ", "", paste(paste("[",paste(sortieLong_sc5_hp$day, sortieLong_sc5_hp$slice)),"]")))
sortieLong_sc5_hp <- mutate(sortieLong_sc5_hp, ds= recode(paste0(day, slice)))

age.l <- c("15-29 yrs", "30-59 yrs", "60-75 yrs")
names(age.l) <- c("1", "2", "3")

sex.l <- c("Male", "Female")
names(sex.l) <- c("1", "2")


p1 <- ggplot(data=sortieLong_sc5_hp %>% filter(facteur=="propHealthy"), aes(x=ds, y=valeur, group=educ, fill=educ, colour=as.character(educ)))
p1 <- p1 + scale_color_discrete(
    name = "Education level",
    labels = c("low", "middle", "high")
  )
p1 <- p1 + theme_bw() + theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))
p1 <- p1 + labs(title="Proportion of healthy by gender, age and education level", y = "Value", x = "Time Slice")
p1 <- p1 + geom_line(size=0.5) + facet_grid( sex ~ age, labeller = labeller(sex = sex.l, age = age.l), scales = "fixed") #show.legend = FALSE
p1 <- p1 + coord_cartesian( ylim = c(0, 0.3))

p2 <- ggplot(data=sortieLong_sc5_hp %>% filter(facteur=="avgOpinion"), aes(x=ds, y=valeur, group=educ, fill=educ, colour=as.character(educ)))
p2 <- p2 + scale_color_discrete(
    name = "Education level",
    labels = c("low", "middle", "high")
  )
p2 <- p2 + theme_bw() + theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))
p2 <- p2 + labs(title="Average Opinion by gender, age and education level", color = "Education level", y = "Value", x = "Time Slice")
p2 <- p2 + geom_line(size=0.5) + facet_grid(sex ~ age, labeller = labeller(sex = sex.l, age = age.l), scales = "fixed") #show.legend = FALSE
p2 <- p2 + coord_cartesian( ylim = c(0.3, 0.7))

propHealtyAndOpinion_higherprop <- plot_grid(
  p1, p2,
  labels = "AUTO", ncol = 1)

print(propHealtyAndOpinion_higherprop)

##ggsave("out/propHealtyAndOpinion_higherprop.pdf")